Efficiency of quantum Monte Carlo impurity solvers for dynamical mean-field theory 
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Since the inception of the dynamical mean-field theory, numerous numerical studies have relied 
on the Hirsch-Fye quantum Monte Carlo (HF-QMC) method for solving the associated impurity 
problem. Recently developed continuous-time algorithms (CT-QMC) avoid the Trotter discretiza- 
tion error and allow for faster configuration updates, which makes them candidates for replacing 
HF-QMC. We demonstrate, however, that a state-of-the-art implementation of HF-QMC (with ex- 
trapolation of discretization At — > 0) is competitive with CT-QMC. A quantitative analysis of 
Trotter errors in HF-QMC estimates and of appropriate At values is included. 

PACS numbers: 71.30.+h, 71.10.Fd, 71.27. +a 



I. INTRODUCTION 

A conventional starting point for the study of strongly 
correlated electron systems is the Hubbard model, which, 
in its single-band version, reads 
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Unfortunately, numerical methods for its direct solution 
are either restricted to onc-dimcnsional cases or suffer, in 
general, from severe finite-size errors and/or sign prob- 
lems. Insight into the physics of higher-dimensional sys- 
tems, thus, requires the use of additional approximations. 
The dynamical mean-field theory (DMFT) neglects inter- 
site correlations by assuming a momentum-independent 
self-energy; it becomes exact in the limit of infinite coor- 
dination number. The DMFT maps the lattice problem 
onto a single-impurity Anderson model (SIAM), supple- 
mented by a self-consistency condition^ Its enormous 
success within the last 15 years would not have been 
possible without the availability of controlled numerical 
solvers for (multi-orbital) SIAMs, in particular, of the 
auxiliary-field Hirsch-Fye quantum Monte Carlo (HF- 
QMC) algorithm^ 

The HF-QMC method discretizes the imaginary-time 
path integral into A time slices of uniform width At = 
(3 /A (at finite temperature T = 1/(3); a Hubbard- 
Stratonovich (HS) transformation replaces the electron- 
electron interaction at each time step by a binary aux- 
iliary field which is sampled by standard Markov Monte 
Carlo (MC) techniques. The numerically costly part in 
each MC step is the update of the Green function (A x A 
matrix), which involves C(A 3 ) operations. HF-QMC is 
numerically exact only in the limit At — > 0; raw results 
(for fixed At) contain the Trotter error, a statistical er- 
ror (which decays as N^ 1 / 2 for N MC sweeps), and - 
in the DMFT context - a convergency error. Keeping 
At fixed for constant accuracy, the computational cost 
increases as T~ 3 , which limits the temperature range ac- 
cessible to HF-QMC. The usage of HF-QMC as a DMFT 
solver requires specialized techniques for Fourier trans- 
forming the time-discretized Green function; violations 



of causality observed in early implementations^ may be 
avoided by stabilizing spline interpolations with analytic 
high-frequency expansions'"^ 

While fundamentally different DMFT solvers such as 
exact diagonalization or renormalization group methods 
often yield very useful information (in particular for the 
one-band case), they are more complementary than gen- 
eral alternatives to HF-QMC, e.g., for obtaining ground 
state or low-frequency results. Thus, much effort is being 
devoted to improving and extending HF-QMC, e.g., by 
incorporating projections to the ground stated and new 
HS decouplings for spin-flip terms in multi-band Hubbard 
models (which otherwise lead to sign problems)' Very re- 
cently, two novel quantum Monte Carlo approaches have 
been formulated for solving the impurity problem in con- 
tinuous imaginary time (CT-QMC), which are based on 
a weak-coupling expansion 8 and an expansion in the im- 
purity hybridization, 9 respectively. Each of these CT- 
QMC methods, which eliminate the Trotter error in HF- 
QMC, is a potential candidate for superseding HF-QMC 
as general-purpose finite-temperature DMFT solver. 

The urgent need for quantitative comparisons between 
the CT-QMC and HF-QMC finite-temperature algo- 
rithms was soon realized. However, a first analysis of 
relative efficiency«ifl suffers from using an inefficient HF- 
QMC-DMFT variant and neglects the essential step of 
extrapolating the HF-QMC results to the limit At = 0. 
It is the purpose of this paper to show that a correction 
of these issues reverses the result: for the given test case 
of a half-filled single-band Hubbard model [Eq. (|TJ)] with 
semi-elliptic density of states at low temperatures and for 
fixed computing resources, extrapolated HF-QMC results 
reach or surpass the precision of the CT-QMC methods. 



II. COMPARISONS AT EQUAL CPU TIME 

In the following, we will extend the performance com- 
parison initiated by Gull et aL— for the case of the 
interaction being chosen equal to the full bandwidth, 
U = W = 4; the scale corresponds to setting t = 1/yZ 
in Eq. (Q]) on the Bcthe lattice in the limit of infinite 
coordination number Z. Unless noted, all CT-QMC re- 
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FIG. 2: (Color online) Comparisons of errors in QMC energy 
estimates: absolute deviations of the data shown in Fig. Q] 
from the reference curve on a logarithmic scale. 



FIG. 1: (Color online) a) Energy versus squared temperature: 
extrapolated HF-QMC results (circles) and estimates (Ref. 
ITol ) from CT-QMC (triangles, squares), b) Data with leading 
T corrections subtracted, plus reference (filled circles). 



suits shown in the following have been extracted from 
Ref. [13 These will be compared with new HF-QMC re- 
sults at equal CPU time: 140 CPU hours on common 
AMD Opteron compute nodes. Specifically, in the case 
of CT-QMC, 20 iterations of 7 hours each have been per- 
formed serially on a single AMD Opteron 244 CPU^ 
In the case of HF-QMC, the same fixed total CPU time 
has been split up among runs for 5-9 different discretiza- 
tions. Keeping, for simplicity, the number of sweeps con- 
stant (10 6 per iteration), this was accomplished by ad- 
justing the number of iterations in a range of 4-8 (which 
appears a bit small for reliable error analysis, but suffi- 
cient for this comparison). The advantage for the HF- 
QMC simulations of being run on a mix of slightly faster 
processors^ 2 - should be almost offset by their overhead of 
parallel execution. 

Figure [Th. compares estimates of the energy, obtained 
from extrapolated HF-QMC, CT-QMC using the hy- 
bridization expansion^ and weak-coupling CT-QMC; 8 
also shown is a solid line which may be considered ex- 
act for the purpose of this comparison. Evidently, all 
QMC methods sustain reasonable accuracy down to fairly 
low temperatures T « IU/200 with moderate computa- 
tional cost. In order to resolve the differences, the lead- 
ing temperature effects have been subtracted in Fig. [TJd. 
At this scale, the CT-QMC results fluctuate visibly in 
a range consistent with their error bars, with possibly 
a slight positive bias for the weak-coupling variant^ 
In contrast, the corresponding HF-QMC results show 
much smaller errors and are in excellent agreement with 
extreme-precision HF-QMC results (filled circles). 

For these reference results, 7 or 8 different discretiza- 
tions using about 30-60 DMFT iterations with up to 
5 • 10 6 Monte Carlo sweeps each have been used for each 
temperature. All raw data points have carefully been 
tested for full convergence; due to the large number of 



iterations, autocorrelation times could be accurately de- 
termined. The HF-QMC runs for the lower temperatures 
have involved matrix sizes of up to A max = 300 with a 
computational cost of about 5000 CPU h per tempera- 
ture; only for T — 1/40 with A max = 320 the effort was 
higher (about 15000 CPU h). The fit curves to this data 
(solid lines in both panels of Fig. [T|) have the form 

E{T) = -0.18051 + 3.941 T 2 - 238 T 4 + 12746 T 6 . 

Here, the leading finite-temperature correction \"]T 2 has 
been obtained from our estimate for the quasiparticle 
weight (at T = 0, U = 4) of Z = 0.2657 ± 0.0006 via 
the relation 7/2 = ir/(3Z); only the other 3 coefficients 
were free fit parameters. For this reason, our ground state 
estimate Eo = —0.18051 ± 0.00003 has about the same 
(small) error bar as the best finite-temperature data. 

A direct comparison of the accuracies of the different 
methods is shown in Fig. [5] The absolute errors (de- 
termined as deviation from the reference curve) of the 
hybridization CT-QMC results (triangles) generally in- 
crease from 10~ 4 at the highest temperatures to 10 -3 at 
lower T; similar behavior (with stronger fluctuations) is 
observed for the weak-coupling CT-QMC data (squares). 
Already the "regular" HF-QMC results (using 140 CPU 
h, open circles) are more accurate by nearly an order 
of magnitude, only to be surpassed by the more costly 
HF-QMC points (solid circles); note that all error esti- 
mates are reliable only down to the precision of the ref- 
erence of about 10~ 5 . In the following, we will highlight 
the methodological ingredients that allow for the high 
demonstrated efficiency of our HF-QMC procedure. 

An essential feature of reliable HF-QMC studies is the 
inclusion of runs with different discretizations At for es- 
timating the impact of the associated systematic error; 
high-precision HF-QMC results can only be obtained by 
extrapolations At — > as illustrated in Fig. [3^: The raw 
HF-QMC results for the energy (circles) vary by about 
0.04 when the discretization is reduced from At = 0.4 
to At = 0.2. As revealed by the extrapolation (dashed 
line), even the best raw results still contain systematic er- 
rors of order 10 -2 ; however, the At error is very regular 
and can essentially be eliminated in a quadratic least- 
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FIG. 3: (Color online) a) HF-QMC estimate of energy at 
T = 1/45 (circles) versus discretization At with extrapolation 
At — > (dashed line), b) Same data after subtraction of 
approximate leading At error versus At 2 in comparison with 
reference HF-QMC data (filled circles and solid line). 

squares fit in At 2 (cf. Sec. IIIip . Already the subtraction 
of the (easily estimated) leading correction (which is lin- 
ear in At 2 ) reduces the errors by one order of magnitude 
as shown in Fig.[3jD; only at this scale, the statistical error 
bars (of about 5 • 10~ 5 ) become visible. Full quadratic 
extrapolation (in At 2 ; dashed line) of the regular HF- 
QMC data results in a total error of 1.5 • 10 -4 (empty 
circle at At 2 = 0), i.e., in a reduction of the initial error 
by two full orders of magnitude. The excellent agreement 
with high-precision HF-QMC data (filled circles and solid 
line) confirms the validity of this procedure. Note that 
already the standard range of discretizations r S [0.2, 0.4] 
in our HF-QMC data implies matrix sizes of up to 226 
(at T = 1/45), vastly greater than typical matrix sizes of 
about 80 for weak-coupling and about 12 for hybridiza- 
tion expansion CT-QMC^ 

An important prerequisite for extrapolation schemes 
such as described above is that all errors in the raw data 
(for fixed At) are well under control. By Fourier trans- 
forming only the difference between the QMC estimated 
Green function and a reference Green function (which is 
exact to order U 2 ), our HF-QMC implementation 3 elim- 
inates At errors beyond the inevitable Trotter contribu- 
tion; at the same time, high-frequency fluctuations are 
reduced. This is illustrated in the comparison of QMC 
self-energy estimates in Fig. |4j Even the HF-QMC result 
for the coarsest discretization At = 0.4 (short-dashed 
line) deviates from the results for At = 0.2 and At — > 
(long-dashed and solid lines, respectively) visibly only for 
^ 3 (magnified in the inset). The extrapolated HF- 
QMC curve, in turn, agrees with the CT-QMC data at 
all frequencies, except for high-frequency fluctuations in 
the hybridization CT-QMC results. 

As a consequence of the minimization of discretization 
errors in our HF-QMC implementation, also fluctuations 
in static quantities such as estimates of the kinetic energy 
are suppressed, as seen in Fig. our HF-QMC results at 
finite discretization (squares/diamonds) follow the refer- 
ence results (dashed lines) with only small scatter; thus, 
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FIG. 4: (Color online) Imaginary part of self-energy on the 
imaginary axis for T = 1/45 as estimated by CT-QMC (Ref. 
10) using the hybridization expansion (triangles) or the weak- 
coupling expansion (circles) and by Hirsch-Fye QMC (dashed 
lines for At = 0.4, At = 0.2, solid line for At -> 0). 



-0.52 



■ 1 1 1 


1 1 


1 1 


"*"- a u. 






O HF-QMC: Ax -> 


***** 




□ HF-QMC: Ax=0. 15 

X HF-QMC (Gull): Ax=0.15 






_ O HF-QMC: At=0.20 




'T 


+ HF-QMC (Gull): Ax=0.20 
I . I 


I 





0.001 0.002 0.003 0.004 

T 2 



FIG. 5: (Color online) Kinetic energy versus squared temper- 
ature (minus leading T corrections): estimates from HF-QMC 
at finite discretization (squares for At = 0.15, diamonds for 
At = 0.2) and after extrapolation At — * (circles); dashed 
and solid lines denote corresponding reference results. Also 
shown: HF-QMC results extracted from Fig. 3 in Ref. 
(crosses); thin lines are guides to the eye only. 

a reliable extrapolation At — > is possible (circles and 
solid line). In contrast, the HF-QMC results from Ref.fTol 
(crosses) show an order of magnitude larger fluctuations 
and do not allow for a precise extrapolation, although 
they have been obtained at higher cost (140 CPU h per 
data point) than each of our finite-AT results (for which 
between 5 and 70 CPU h were used, depending on tem- 
perature}. Evidently, the HF-QMC implementation used 
in Ref. ll(J is much less efficient. 

As seen in Fig. [6l the estimates of the kinetic en- 
ergy obtained from weak-coupling CT-QMC and At- 
extrapolated regular HF-QMC agree within their (sim- 
ilar) error bars with the reference HF-QMC results; the 
extrapolation of the latter for T — > (solid line) is also in 
excellent agreement with low-T hybridization CT-QMC 
data of Ref. [9| (double-dashed line). The regular hy- 
bridization CT-QMC data (triangles) shows by a factor 
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FIG. 6: (Color online) Kinetic energy versus squared tem- 
perature: estimates at fixed CPU time from CT-QMC (Ref. 
Ilpf ) (triangles, squares) and HF-QMC (empty circles) after 
subtracting leading T corrections. Also shown are higher- 
precision results from HF-QMC (filled circles and solid line) 
and hybridization CT-QMC (double-dashed line). 



of about 2 larger deviations from the reference (this devi- 
ation is systematic and far beyond the stated error bars). 
So, all 3 QMC methods show roughly the same efficiency 
in estimating this observable. 



III. QUANTITATIVE STUDY OF HF-QMC 
TROTTER ERRORS 

As we have established above for the specific test case 
defined in Ref. 10, HF-QMC results should be extrapo- 
lated to At = 0, if possible. Many HF-QMC studies, 
however, include simulations only for a single value of 
At. The reliability of such studies depends on whether 
the chosen values of At can be regarded as "small" . Up 
to now, many different criteria have been used in the lit- 
erature, such as (for t* = 1) At < 0.25 or At < 1/(5Z7) 
etc^ We have invested significant computing resources 
for determining the leading and subleading Trotter errors 
in HF-QMC estimates of selected observables throughout 
the (paramagnetic) phase space of model fl} . Figure [7| 
shows the coefficients in an expansion of the form 

E{U,T;At) « E < - a \U,T)+E^(U,T)AT 2 +E^(U,T)^T i 

versus interaction U for a range of temperatures extend- 
ing from a value slightly above the critical Mott end point 
(U w 4.67, T sa 0.055) to much smaller values. For 
weak interaction U < 2, both the leading correction E^ 2 > 
(lower panel) and the subleading term E^ (upper panel) 
scale asymptotically as U 2 (double-dashed lines); thus, 
in this regime the Trotter error is indeed proportional to 
(C/At) 2 . For larger interactions, however, the error grows 
much slower; it even decreases as a function of U in the 
strongly correlated metallic phase (4 < U < 5). Only in 
this range, a temperature dependence of the coefficients 
becomes clearly visible. Finally, in the insulating phase 
(U > 5), the leading coefficient E^ grows essentially 
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FIG. 7: (Color online) Coefficients of Trotter errors in HF- 
QMC estimates of the energy (see text). 




FIG. 8: (Color online) Coefficients of Trotter errors in HF- 
QMC estimates of kinetic energy. Inset: relative coefficients. 

linearly with U, so that AtVU would have to be kept 
constant for constant Trotter error. 

The corresponding coefficients in the HF-QMC Trot- 
ter error expansion for the kinetic energy show markedly 
different behaviors in Fig. [5] Here, the leading coefficient 
-^kin sca l es as U 3 for U < 3; E^ is too noisy for this kind 
of analysis. Both coefficients are strongly enhanced near 
the Mott transitions, with a pole-like structure in 
Again, the temperature dependences are appreciable only 
in this range. Somewhat surprisingly, the Trotter errors 
decay quickly in the insulating phase, even the relative 
errors (see inset of Fig. [5]); this observation may be spe- 
cific to our implementation. 

In the case of the double occupancy D — {n^riii) (see 
Fig. [5]), the Trotter error starts off as At 2 £7 for very 
small interactions U < 1. It then decays with a sign 
change just below the Mott transition. In the insulating 
phase, both coefficients are nearly constant; a decay of 
for extremely large U is visible in the inset. 
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FIG. 9: (Color online) Coefficients of Trotter errors in HF- 
QMC estimates of double occupancy. 



IV. CONCLUSION 

We have compared published results 1 ^ of two different 
continuous-time QMC methods with data obtained from 
a Hirsch-Fye QMC implementation by extrapolating the 
discretization At — > 0. HF-QMC involved typical matrix 
sizes more than twice as large as weak-coupling CT-QMC 
and nearly 20 times larger than hybridization CT-QMC 
(cf. Fig. 2 inRef.[l0|). If this was the only relevant factor, 
HF-QMC should (due to the cubic scaling of the runtime 
of each method with the matrix size) be less efficient 
by one order of magnitude than the weak-coupling and 
by four orders of magnitude than the hybridization ex- 
pansion variant. However, at (roughly) equal numerical 
effort, the HF-QMC results are about equally accurate 
for the kinetic energy and an order of magnitude more 
accurate for the energy; the latter translates^ into two 
orders of magnitude higher efficiency of HF-QMC. Ap- 
parently, much larger fluctuations in the CT-QMC meth- 
ods outweigh the advantage of using smaller matrices (in 



the temperature range included in this comparison). An 
important factor is certainly the computation of Green 
functions which is (on the discretization grid) trivial for 
HF-QMC (with A measurements for each time difference 
Ti — Tj resulting from each sweep) while it requires sig- 
nificant additional efforts in both CT-QMC methods^ 

We have also quantified contributions to the Trotter er- 
ror in HF-QMC estimates of energetics which show quite 
complex behavior. A rule for "small enough" values of 
At which holds approximately for large U can be ob- 
tained by requiring the quartic corrections (for E or D) to 
be smaller than the quadratic ones: At -C l/\/0.28{7f*. 
The observed small temperature dependence of the coef- 
ficients could be used for cheap high-precision studies at 
low T (by extrapolating the coefficients and performing 
low-T simulations only for large At). A more surpris- 
ing insight is that the Trotter errors in D are unusually 
small for U ~ 4.8 and T < 40; the same holds true 
for E^ so that HF-QMC is particularly precise in this 
range. More generally, HF-QMC remains the method of 
choice for obtaining precise energetics, even extrapolated 
to the ground stated 

The CT-QMC methods are certainly important and 
promising additions to the portfolio of DMFT impurity 
solvers. They allow for direct simulations at extremely 
low temperatures without requiring extrapolations in At 
and/or T and appear to offer greater flexibility in the 
choice of the Hamiltonian without incurring significant 
sign problems. However, as shown in this study, HF- 
QMC (with extrapolation At — * 0) can be equally or 
even more efficient and should not be discarded (yet). 
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